module Population
!  File: Popuplation.f90
!   Purpose: Contains the poplation models
!      Models
!         1. Separable
!         2. VPA
!         3. ISVPA
!         4. PS (Pope & Shepherd)
!   Note:
!      Model Input
!         Observed Catch in numbers (canum)
!         Parameters
!            F Year factor
!            Selection
!            Natural mortality
!            N Last year
!            N oldest age
!   History: When        Who        What
!            25-04-07    Skalevik   Created
   use OptimisationTypes
   use ObsData
   use ObsModel
   use Parameter
   use Equations
   use matrix
   use HTypes
   use Index
   implicit none

!logical, parameter :: g_UseISVPAScaling = .true.
!logical, parameter :: g_UseNormalising = .false.

contains


subroutine Population_Allocate(Opt)
!   Purpose: Initializes the Fay and Nay matrixes
!   History: When        Who        What
!            25-04-07    Skalevik   Created
!            08-04-08    Skalevik   Enhanced - allocates space for next year 
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iAge, iYear, NYears, NAges
   ! Reallocate population
   allocate(Opt%PopModel)
   NYears = NumYears(Opt%IndexInput) + 1
   NAges = NumAges(Opt%IndexInput)
   if ((0 < NYears) .and. (0 < NAges)) then
      allocate(Opt%PopModel%Nay(0 : NYears - 1, 0 : NAges - 1))
      call InitMatrix(Opt%PopModel%Nay, 0.0d0)
      allocate(Opt%PopModel%Fay(0 : NYears - 1, 0 : NAges - 1))
      call InitMatrix(Opt%PopModel%Fay, 0.0d0)
      allocate(Opt%PopModel%May(0 : NYears - 1, 0 : NAges - 1))
      call InitMatrix(Opt%PopModel%May, 0.0d0)
      allocate(Opt%PopModel%InternalCatchISVPA(0 : NYears - 1, 0 : NAges - 1))
      call InitMatrix(Opt%PopModel%InternalCatchISVPA, 0.0d0)
      allocate(Opt%PopModel%StockSummary(0 : NYears - 1, ssFirst:ssLast))
      call InitMatrix(Opt%PopModel%StockSummary, 0.0d0)
   end if
end subroutine Population_Allocate


subroutine DeallocatePopulation(Opt)
!   Purpose: Deallocates the Fay and Nay matrixes
!   History: When        Who        What
!            02-06-07    Skalevik   Created

   implicit none
   type(TOptimisation), intent(inout) :: Opt
   if (associated(Opt%PopModel)) then
      if (associated(Opt%PopModel%Nay)) then
         deallocate(Opt%PopModel%Nay)
      end if
      if (associated(Opt%PopModel%Fay)) then
         deallocate(Opt%PopModel%Fay)
      end if
      if (associated(Opt%PopModel%May)) then
         deallocate(Opt%PopModel%May)
      end if
      if (associated(Opt%PopModel%InternalCatchISVPA)) then
         deallocate(Opt%PopModel%InternalCatchISVPA)
      end if
      if (associated(Opt%PopModel%StockSummary)) then
         deallocate(Opt%PopModel%StockSummary)
      end if
      deallocate(Opt%PopModel)
   end if
end subroutine DeallocatePopulation


subroutine NormalizeSelection(Opt, ScaleYFactor)
!   Purpose: Selection and annual F share a scaling factor. The selection is scaled 
!            so that the average over a reference age range is one. That gives values 
!            of the annual F scaled as the normally reported average F.
!   Note:
!   History: When        Who        What
!            14-06-07    Skalevik   Created
   implicit none
   type(TOptimisation), intent(in) :: Opt
   logical, intent(in) :: ScaleYFactor
   double precision :: SumRefSel, ScaleFac, NewYFactor
   integer :: iYear, iAge, iFirstRef, iLastRef, NumRefAges
   iFirstRef = Opt%IndexInput%FirstRefAge - Opt%IndexInput%FirstAge
   iLastRef = Opt%IndexInput%LastRefAge - Opt%IndexInput%FirstAge
   do iYear = 0, NumYears(Opt%IndexInput) - 1
      SumRefSel = 0
      ScaleFac = 0
      do iAge = iFirstRef, iLastRef
         SumRefSel = SumRefSel + PValue(Opt, cSelection, 0, iYear, iAge)
      end do
      NumRefAges = iLastRef - iFirstRef + 1
      if (0 /= SumRefSel) then
         ScaleFac = NumRefAges / SumRefSel
         ! Scale total age-range selection
         do iAge = 0, NumAges(Opt%IndexInput) - 1
            call SetPValue(Opt, cSelection, 0, iYear, iAge, PValue(Opt, cSelection, 0, iYear, iAge) * ScaleFac)
         end do
         if(ScaleYFactor) then
            NewYFactor = PValue(Opt, cFYearFactor, 0, iYear, 0) / ScaleFac
            call SetPValue(Opt, cFYearFactor, 0, iYear, 0, NewYFactor)
         end if
      end if  
   end do
end subroutine NormalizeSelection


function FyScaleForISVPA(Opt)
!   Purpose: Calculate the scale factor for F yearfactor into isvpa routine
!   History: When        Who        What
!            25-08-08    Skalevik   Created
   implicit none
   type(TOptimisation), intent(in) :: Opt
   double precision :: FyScaleForISVPA
   integer :: iAge, iYear, iFirstRef, iLastRef
   double precision :: SumRefSel, SumAllSel
   
   FyScaleForISVPA = 1.0d0
   iFirstRef = Opt%IndexInput%FirstRefAge - Opt%IndexInput%FirstAge
   iLastRef = Opt%IndexInput%LastRefAge - Opt%IndexInput%FirstAge
   SumRefSel = 0.0d0
   SumAllSel = 0.0d0
   iYear = NumCatchYears(Opt%IndexInput) - 1
   do iAge = 0, NumTrueAges(Opt%IndexInput) - 1
      if ((iAge >= iFirstRef) .and. (iAge <= iLastRef)) then
         SumRefSel = SumRefSel + PValue(Opt, cSelection, 0, iYear, iAge)
      endif
      SumAllSel = SumAllSel + PValue(Opt, cSelection, 0, iYear, iAge)
   end do
   if ((SumAllSel /= 0.0d0) .and. (SumRefSel /= 0.0d0)) then
      FyScaleForISVPA = SumAllSel / SumRefSel 
   endif
end function FyScaleForISVPA


subroutine SetFay(Opt, iYear, iAge, FVal)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iYear, iAge
   double precision :: FVal
!   if(FVal < 0) then
!      if (0.000001 < abs(FVal)) then
!         call MessageOut("Error: Negative F value")
!      end if
!      FVal = 0
!   endif
   Opt%PopModel%Fay(iYear, iAge) = FVal
   if ((FVal < 0.0d0) .and. (FVal /= -1.0d0)) then
      write(*,*) "FATAL ERROR in SetFay (I) F < 0"
      call MessageOut(Opt, " F="  // floattostr(FVal))
      stop
   end if
end subroutine SetFay


subroutine SetNay(Opt, iYear, iAge, NVal)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iYear, iAge
   double precision :: NVal
   !if(NVal < 0) then
   !   write(*,*) "setnay error"
   !   write(*,*)  iYear, iage, NVal 
   !   stop
   !endif
   Opt%PopModel%Nay(iYear, iAge) = NVal
end subroutine SetNay


subroutine CalculatePopulation(Opt, Finished)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   logical, intent(in) :: Finished
   integer :: iYear, iAge

   !call NormalizeSelection(Opt, .false.)
   if (NumYears(Opt%IndexInput) < 1) then
      call MessageOut(Opt, "# Years must be greater than 0.")
      return
   endif
   if (NumTrueAges(Opt%IndexInput) < 1) then
      call MessageOut(Opt, "# True ages must be greater than 0.")
      return
   endif
   select case(Opt%IndexInput%ModelType)
      case(pmtSep)
         call CalculateSeparablePopulation(Opt)
      case(pmtVPA) 
         call CalculateVPAPopulation(Opt)
      case(pmtISVPA) 
         call CalculateISVPAPopulation(Opt, Finished)
   end select
   call OnCalculateFinished(Opt)
end subroutine CalculatePopulation


subroutine NormalizeSelectionLocally(Opt, NewSel)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   double precision, pointer :: NewSel(:, :)
   integer :: iYear, iAge, iFirstRef, iLastRef, NumRefAges, Num_TrueAges, Num_Years
   double precision :: SumRefSel, ScaleFac, Sel
   Num_TrueAges = NumTrueAges(Opt%IndexInput)
   Num_Years = NumYears(Opt%IndexInput)
   ! Normalize selection locally
   iFirstRef = Opt%IndexInput%FirstRefAge - Opt%IndexInput%FirstAge
   iLastRef = Opt%IndexInput%LastRefAge - Opt%IndexInput%FirstAge
   do iYear = 0, Num_Years - 1
      SumRefSel = 0
      ScaleFac = 0
      do iAge = iFirstRef, iLastRef
         Sel = PValue(Opt, cSelection, 0, iYear, iAge)
         SumRefSel = SumRefSel + Sel
      end do
      NumRefAges = iLastRef - iFirstRef + 1
      if (0 /= SumRefSel) then
         ScaleFac = NumRefAges / SumRefSel
         ! Scale total age-range selection
         do iAge = 0, Num_TrueAges - 1
            Sel = PValue(Opt, cSelection, 0, iYear, iAge)
            Sel = Sel * ScaleFac
            NewSel(iYear, iAge) = Sel
         end do
      end if  
   end do
end subroutine NormalizeSelectionLocally

 
!function GetPopulationFromFleet(Opt, N0, iYear, iAge)
!   implicit none
!   type(TOptimisation), intent(inout) :: Opt
!   double precision, intent(out) :: N0
!   integer, intent(in) :: iYear, iAge
!   logical :: GetPopulationFromFleet
!   integer :: NumDataSources, NValidIndexes, iDataSource, NValidFleets
!   double precision :: ObsVal, E_ObsVal, q, E_N0, Z, M, F
!	type(TDataSource), pointer :: pObsData
!   NumDataSources = GetNumObservedDataElms(Opt, odFleet) 
!   NValidFleets = 0
!   E_N0 = 0
!   N0 = 0;
!   M = PValue(Opt, cNatMort, 0, iYear, iAge) 
!   F = 0.0d0 ! F assumed to be zero
!   if (0 < NumDataSources) then
!      do iDataSource = 0, NumDataSources - 1
!        pObsData => GetDataSource(Opt, odFleet, iDataSource)
!         ObsVal = GetObsDataValue(Opt, odFleet, iDataSource, iYear, iAge)
!   	   q = PValue(Opt, cCBL_Age, iDataSource, iYear, iAge) 
!         if ((0.0 < ObsVal) .and. (0 < q)) then
!            NValidFleets = NValidFleets + 1
!            E_ObsVal = E_ObsVal + ObsVal
!            Z = F * pObsData%FProp + M * pObsData%MProp 
!           E_N0 = E_N0 + ObsVal / q * SafeExp(Z)
!         endif
!      enddo 
!   endif
!   GetPopulationFromFleet = 0 < NValidFleets
!   if (GetPopulationFromFleet) then
!      N0 = E_N0 / NValidFleets
!   endif
!end function GetPopulationFromFleet


subroutine EstimateVPAInOldestTrueAge(Opt, iYear, iYearLastValidF, LastValidFInVPAOldestAge)
   type(TOptimisation), intent(inout) :: Opt
   integer, intent(in) :: iYear, iYearLastValidF
   double precision, intent(inout) :: LastValidFInVPAOldestAge
   double precision :: E_F_Ay0, E_F_Ayn, F_Ay0, F_Ayn, Fay0, Fayn, FRate
   integer :: iAge, Nay, Num_TrueAges, Num_Years, iTrueAge, iFirstRef, iLastRef
   
   Num_TrueAges = NumTrueAges(Opt%IndexInput)
   Num_Years = NumYears(Opt%IndexInput)
   E_F_Ay0 = 0.0
   E_F_Ayn = 0.0
   Nay = 0
   iFirstRef = Opt%IndexInput%FirstRefAge - Opt%IndexInput%FirstAge
   iLastRef = Opt%IndexInput%LastRefAge - Opt%IndexInput%FirstAge
   !iTrueAge = Num_TrueAges - 1
   if (iFirstRef < iLastRef) then
      do iAge = iFirstRef, iLastRef - 1
         Fay0 = Opt%PopModel%Fay(iYear, iAge)      
         Fayn = Opt%PopModel%Fay(iYearLastValidF, iAge)      
         if ((Fay0 <= 0.0) .or. (Fayn <= 0.0)) cycle
         E_F_Ay0 = E_F_Ay0 + Fay0
         E_F_Ayn = E_F_Ayn + Fayn
         Nay = Nay + 1
      enddo
   endif
   FRate = 1.0d0
   if ((E_F_Ay0 > 0.0d0) .and. (E_F_Ayn > 0.0d0)) then
      FRate = ((E_F_Ay0 / Nay) / (E_F_Ayn / Nay))
   endif
   
   ! The F_ayn estimation
   ! F = F_ayn * Rate     
   F_Ayn = LastValidFInVPAOldestAge;
   LastValidFInVPAOldestAge = F_Ayn * FRate
end subroutine EstimateVPAInOldestTrueAge


subroutine CalculateSeparablePopulation(Opt)
!   Purpose: Calculates the separable population,
!            Fishing mortalities is calculated by the product of selection Say and yearfactor Fy
!   Note:
!   History: When        Who        What
!            25-04-07    Skalevik   Created
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iYear, iAge, iCohort, NumCohorts, &
      Num_TrueAges, Num_Years, iYearEnd, iAgeEnd
   double precision :: F, S, Fy, M, Z, N0, N1
   double precision, pointer :: NewSel(:, :)

   Num_TrueAges = NumTrueAges(Opt%IndexInput)
   Num_Years = NumYears(Opt%IndexInput)

   allocate(NewSel(0 : Num_Years - 1, 0 : Num_TrueAges - 1))
   call NormalizeSelectionLocally(Opt, NewSel)
   
   ! Trace the cohorts backwards
   NumCohorts = Num_TrueAges + Num_Years - 1
   do iCohort = 0, NumCohorts - 1
      ! Find the indexes at the end of the cohort, starting with the youngest one.
      if (iCohort <= Num_TrueAges - 1) then
         iYearEnd = Num_Years - 1
         iAgeEnd = iCohort
      else
         iAgeEnd = Num_TrueAges - 1
         iYearEnd = Num_Years - 1 - (iCohort - Num_TrueAges + 1)
      endif
      iYear = iYearEnd
      iAge = iAgeEnd
      do 
         ! Pick N1:
         if (iYear == Num_Years - 1) then
            N1 = PValue(Opt, cNValues, 0, Num_Years - 1, iAge)
         else if (iAge == Num_TrueAges - 1) then
            N1 = PValue(Opt, cNValues, 0, iYear, Num_TrueAges - 1)
         else
            N1 = Opt%PopModel%Nay(iYear + 1, iAge + 1)
         endif
         ! Compute F and Z:
         M = PValue(Opt, cNatMort, 0, iYear, iAge) 
         Fy = PValue(Opt, cFYearFactor, 0, iYear, 0)
         S = NewSel(iYear, iAge)
         F = Fy * S
         Z = F + M 
         ! Compute N0:
         if (Z < 80) then
            N0 = N1 * Exp(Z)
         else
            N0 = 0
         endif
         ! Set F and N0:
         call SetFay(Opt, iYear, iAge, F)
         if(N0 < 0.0d0) then
            write(*,*) "fatal error, n0=", N0
            call system("pause")
         endif 
         call SetNay(Opt, iYear, iAge, N0)
         if ((iYear == 0) .or. (iAge == 0)) exit ! Cohort finished
         ! Prepare for the next element in the cohort
         iYear = iYear - 1
         iAge = iAge - 1
      enddo
   enddo
   deallocate(NewSel)
end subroutine CalculateSeparablePopulation


subroutine CalculateVPAPopulation(Opt)
!   Purpose: Calculates the VPA population by tuning survivors.
!            Fishing mortalities is calculated backward over years
!            The stock number is approximately resolved by Popes formula back in time
!   Note:
!   History: When        Who        What
!            27-04-07    Skalevik   Created
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iYear, iAge, iYearEnd, iAgeEnd, iYearForward, iAgeForward, Num_TrueAges, &
      Num_Years, iCohort, NumCohorts, iYearLastValidF, EndFlag;
   double precision :: C, N0, N1, F, M, Z, N0byN1, LastValidFInVPAOldestAge
   logical :: InLastElm, InOldestAge, InLastYear
   character(len = 256) :: str

   Num_TrueAges = NumTrueAges(Opt%IndexInput)
   Num_Years = NumCatchYears(Opt%IndexInput)

   if (Num_Years == 0) return
   if (Num_TrueAges == 0) return
   
   LastValidFInVPAOldestAge = -1
   iYearLastValidF = Num_Years - 1
   NumCohorts = Num_TrueAges + Num_Years - 1
   
   ! Calculate F_VPA and N_VPA
   do iCohort = 0, NumCohorts - 1

      ! Find the indexes at the end of the cohort, starting with the youngest one.
      if (iCohort <= Num_TrueAges - 1) then
         ! In last year
         iYearEnd = Num_Years - 1
         iAgeEnd = iCohort
      else
         ! In oldest age
         iAgeEnd = Num_TrueAges - 1
         iYearEnd = Num_Years - 1 - (iCohort - Num_TrueAges + 1)
      endif
      EndFlag = PFlag(Opt, cNValues, 0, iYearEnd, iAgeEnd)

      iYear = iYearEnd
      iAge = iAgeEnd
      N1 = 0
      
      do 
         ! Iterate each element in cohort backwards
         ! N1 > 0 means that the chohort has been started (initialized N0 in the previous iteration)
         InOldestAge = iAge == Num_TrueAges - 1
         InLastYear = iYear == Num_Years - 1
         InLastElm = InLastYear .or. InOldestAge

         C = max(0.0d0, GetObsDataValue(Opt, odCaNum, 0, iYear, iAge))
         M = PValue(Opt, cNatMort, 0, iYear, iAge)
         
         if (InOldestAge .and. (.not. InLastYear) .and. (EndFlag == flSpecialCode)) then
            ! N0(F) calculation - Calculate N0 from estimated F0
            if (0.0d0 < LastValidFInVPAOldestAge) then
               call EstimateVPAInOldestTrueAge(Opt, iYear, iYearLastValidF, LastValidFInVPAOldestAge)
            endif
            F = LastValidFInVPAOldestAge
            if (F <= 0.0d0) then
               N0 = 0
               !write(*,*) "LastValidFInVPAOldestAge: ", LastValidFInVPAOldestAge
               !call system("pause")
            else
               Z = F + M 
               N0 = C * Z / F / (1 - Exp(-Z))    
            endif
         else
            ! Ordinary VPA. N0(N1). Calculate N0 from N1
            if(InLastElm) then
               N1 = PValue(Opt, cNValues, 0, iYear, iAge)
            else
               N1 = Opt%PopModel%Nay(iYear + 1, iAge + 1)
            endif
            ! Popes:
            N0 = N1 * Exp(M) + C * Exp(M * 0.5)
            ! Fishing mortality from N0
            if((N1 < 0.0d0) .or. (N0 < 0.0d0)) then
               write(*,*) "D Error M, C, N0, N1: ", M, C, N0, N1
            endif
            if ((N0 == 0) .or. (N1 == 0)) then
               F = -1
            else  
               N0byN1 = N0 / N1
               F = log(N0byN1) - M
               ! Precision correction:
               if (F < 0.0d0) then
                  F = 0.0d0
               endif
            endif
            if (InLastYear .and. (F > 0.0d0)) then
               LastValidFInVPAOldestAge = F
            endif
         endif
         if (InOldestAge .and. (F > 0.0)) then
            iYearLastValidF = iYear
         endif

         call SetNay(Opt, iYear, iAge, N0)
         call SetFay(Opt, iYear, iAge, F)

         if ((iYear == 0) .or. (iAge == 0)) exit ! Cohort finished
         ! Prepare for the next element in the cohort
         iYear = iYear - 1
         iAge = iAge - 1
      enddo
   enddo
end subroutine CalculateVPAPopulation


subroutine CalculateISVPAPopulation(Opt, Finished)
!   Purpose: Calculates the ISVPA population
!   Note:    Formula summary in this procedure:
!            Separabel fishing mortality:    F = Fy * S
!            Fishing mort. derived from N    F = LN(N0 / N1) - M
!                                            N0 = N1 * Exp(Z) 
!            Catch eq.                       I)  N0 = C * Z / F / (1 - exp(-Z)) 
!                                            II) C = F * N0 * (1 - exp(-Z)) / Z
!            Popes eq.                       N0 = N1 * Exp(M) + C * Exp(M / 2)  
!   Note:
!   History: When        Who        What
!            27-04-07    Skalevik   Created
!            15-06-07    Skalevik   Enhanced - Leave the last year factor unchanged.
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   logical, intent(in) :: Finished
   double precision :: Fy_New, S_New, Fy_Old, S_Old, F_VPA, N_VPA, F_SEP, M, C_Old, C_New, N0, N1, Z, &
      Res, ResSqr, TotResSqr, F, Fy, S
   double precision, pointer :: Fy_Source(:) 
   double precision, pointer :: S_Source(:, :)  ! Year group - age
   double precision, pointer :: Fy_Smooth(:) 
   double precision, pointer :: S_Smooth(:, :)  ! Year group - age
   double precision, pointer :: ResArr(:, :) 
   double precision, pointer :: ResSqrArr(:, :) 
   integer, pointer :: YearGroups(:)
   integer :: Flag, NumCohorts, NSplits, iCohort, iAge, iYear, iYearEnd, iAgeEnd, Num_TrueAges, &
      Num_Years, NElms, iSplit, iYearGroupStart, iYearGroupEnd
   double precision :: TempSum, TempValue, TempAverage, Convergence, FyScaled
   integer :: TempNum
   logical :: EndOfCohort
   double precision :: LastRes(0 : 3 - 1) 
   type(TParamElm), pointer :: pParamElm

   double precision, parameter :: ISVPA_Limit = 1e-5

   Num_TrueAges = NumTrueAges(Opt%IndexInput)
   Num_Years = NumCatchYears(Opt%IndexInput)

   if (Num_Years == 0) return
   if (Num_TrueAges == 0) return
   if (Finished) then
      allocate(ResArr(0 : Num_Years - 1, 0 : Num_TrueAges - 1)) 
      allocate(ResSqrArr(0 : Num_Years - 1, 0 : Num_TrueAges - 1)) 
   endif
   NumCohorts = Num_TrueAges + Num_Years - 1

   ! Find year groups by splits in the selection 
   NSplits = 1
   if (Num_Years > 0) then
      do iYear = 1, Num_Years - 1
         Flag = PFlag(Opt, cSelection, 0, iYear, 0)
         if (Flag == flEstimate) then
            NSplits = NSplits + 1
         endif
      enddo
   endif
   allocate(YearGroups(0 : NSplits - 1))
   YearGroups(0) = 0
   iSplit = 1
   if (Num_Years > 0) then
      do iYear = 1, Num_Years - 1
         Flag = PFlag(Opt, cSelection, 0, iYear, 0)
         if (Flag == flEstimate) then
            YearGroups(iSplit) = iYear
            iSplit = iSplit + 1
         endif
      enddo
   endif
   
   allocate(S_Source(0 : NSplits - 1, 0 : Num_TrueAges - 1))
   allocate(Fy_Source(0 : Num_Years - 1))
   
   ! Initialize S_Source 
   do iSplit = 0, size(YearGroups) - 1
      do iAge = 0, Num_TrueAges - 1
         S_Source(iSplit, iAge) = 1.0d0 / Num_TrueAges
         !S_Source(iSplit, iAge) = PValue(Opt, cSelection, 0, YearGroups(iSplit), iAge)
      enddo
   enddo
   ! Initialize Fy_Source 
   FyScaled = PValue(Opt, cFYearFactor, 0, Num_Years - 1, 0) * FyScaleForISVPA(Opt)
   do iYear = 0, Num_Years - 1
      Fy_Source(iYear) = FyScaled
      !Fy_Source(iYear) = PValue(Opt, cFYearFactor, 0, iYear, 0) 
   enddo
   
   if(Finished) then
      do iYear = 0, Num_Years - 1
         do iAge = 0, Num_TrueAges - 1
            ResArr(iYear, iAge) = 0.0d0
            ResSqrArr(iYear, iAge) = 0.0d0
         enddo
      enddo
   endif
   
   allocate(S_Smooth(0 : NSplits - 1, 0 : Num_TrueAges - 1))
   allocate(Fy_Smooth(0 : Num_Years - 1))

   LastRes(0 : 3 -1) = 0
   do
      ! Calculate F_VPA and N_VPA
      do iCohort = 0, NumCohorts - 1
         ! Find the indexes at the end of the cohort, starting with the youngest one.
         if (iCohort <= Num_TrueAges - 1) then
            iYearEnd = Num_Years - 1
            iAgeEnd = iCohort
         else
            iAgeEnd = Num_TrueAges - 1
            iYearEnd = Num_Years - 1 - (iCohort - Num_TrueAges + 1)
         endif
         iYear = iYearEnd
         iAge = iAgeEnd
         do 
            do iSplit = size(YearGroups) - 1, 0, -1 
               if (iYear >= YearGroups(iSPlit)) exit
            enddo
            Fy_Old = Fy_Source(iYear)
            S_Old = S_Source(iSplit, iAge)
            !Fy_Old = PValue(Opt, cFYearFactor, 0, iYear, 0) 
            !S_Old = PValue(Opt, cSelection, 0, iYear, iAge) 
            M = PValue(Opt, cNatMort, 0, iYear, iAge) 
            C_Old = GetObsDataValue(Opt, odCaNum, 0, iYear, iAge)
            ! Default mortality is separable: 
            ! F = Fy * S
            F_VPA = Fy_Old * S_Old
            N0 = 0.0d0
            if ((F_VPA > 0) .and. (C_Old >= 0)) then
               Z = F_VPA + M 
               ! Default N is derived from catch by catch eq:
               ! N = C * Z / F / (1 - exp(-Z))
               N0 = C_Old * Z / F_VPA / (1 - Exp(-Z))    
               EndOfCohort = (iYear == iYearEnd) .or. (iAge == iAgeEnd)
               if (.not. EndOfCohort) then
                  N1 = Opt%PopModel%Nay(iYear + 1, iAge + 1)
                  if (N1 > 0) then
                     ! If not end of cohort AND N1 > 0 derive N from catch (popes):
                     ! N0 = N1 * Exp(M) + C * Exp(M / 2)
                     N0 = N1 * Exp(M) + C_Old * Exp(M * 0.5) 
                  endif 
                  if ((C_Old > 0) .and. (N0 > 0) .and. (N1 > 0)) then
                     ! If not end of cohort AND Catch > 0 AND N0 > 0 AND N1 > 0, derive F from N:
                     ! F = LN(N0 / N1) - M
                     F_VPA = Log(N0 / N1) - M
                     if (F_VPA < 0.0d0) then
                        F_VPA = 0.0d0
                     endif
                  endif 
               endif
            endif
            call SetFay(Opt, iYear, iAge, F_VPA)
            call SetNay(Opt, iYear, iAge, N0)
            
            if ((iYear == 0) .or. (iAge == 0)) exit ! Cohort finished
            ! Prepare for the next element in the cohort
            iYear = iYear - 1
            iAge = iAge - 1
         enddo
      enddo
      
      ! Calculated splitted selection patterns
      do iAge = 0, Num_TrueAges - 1
         do iSplit = 0, NSplits - 1
            if ((0 < iAge) .and. (iAge == Num_TrueAges - 1)) then
               TempAverage = S_Smooth(iSplit, iAge - 1)
            else
               iYearGroupStart = YearGroups(iSplit)
               if (iSplit == (NSplits - 1)) then
                  iYearGroupEnd = Num_Years - 1
               else
                  iYearGroupEnd = YearGroups(iSplit + 1) - 1
               endif
               TempSum = 0
               TempNum = 0
               do iYear = iYearGroupStart, iYearGroupEnd
                  F_VPA = Opt%PopModel%Fay(iYear, iAge)
                  !Fy_Old = PValue(Opt, cFYearFactor, 0, iYear, 0) 
                  Fy_Old = Fy_Source(iYear)
                  if (Fy_Old /= 0) then
                     TempValue = F_VPA / Fy_Old
                     TempSum = TempSum + TempValue
                     TempNum = TempNum + 1
                  endif
               enddo
               if (TempNum > 0) then
                  TempAverage = TempSum / TempNum
               else 
                  TempAverage = 0
               endif
            endif
            S_Smooth(iSplit, iAge) = TempAverage
         enddo
      enddo
      
      ! Calculate new S
      do iSplit = 0, NSplits - 1
         ! Calculate sum
         TempSum = 0
         do iAge = 0, Num_TrueAges - 1
            TempSum = TempSum + S_Smooth(iSplit, iAge)
         enddo
         ! Normalize to sum
         if (TempSum > 0) then
            do iAge = 0, Num_TrueAges - 1
               S_Smooth(iSplit, iAge) = S_Smooth(iSplit, iAge) / TempSum
            enddo
         endif
      enddo
      
      ! Calculate new Fy (except for the last year = free parameter)
      do iYear = 0, Num_Years - 1 - 1
         do iSplit = size(YearGroups) - 1, 0, -1 
            if (iYear >= YearGroups(iSPlit)) exit
         enddo
         TempSum = 0
         TempNum = 0
         do iAge = 0, Num_TrueAges - 1
            F_VPA = Opt%PopModel%Fay(iYear, iAge)
            S_New = S_Smooth(iSplit, iAge)
            if (0 < S_New) then
               TempValue = F_VPA / S_New 
               TempSum = TempSum + TempValue
               TempNum = TempNum + 1
            endif
         enddo
         if (TempNum > 0) then
            TempAverage = TempSum / TempNum 
         else
            TempAverage = 0
         endif
         ! Smooth Fy
         Fy_Smooth(iYear) = TempAverage
         !Fy_Old = PValue(Opt, cFYearFactor, 0, iYear, 0) 
         Fy_Old = Fy_Source(iYear)
         Fy_Smooth(iYear) = (1 - Fy_Smooth_Factor) * Fy_Old + Fy_Smooth_Factor * Fy_Smooth(iYear)
      enddo
      Fy_Smooth(Num_Years - 1) = Fy_Source(Num_Years - 1) !PValue(Opt, cFYearFactor, 0, Num_Years - 1, 0) 

      ! Smooth S
      do iSplit = 0, NSplits - 1
         ! Calculate sum
         !iYear = YearGroups(iSplit)
         do iAge = 0, Num_TrueAges - 1
            !S_Old = PValue(Opt, cSelection, 0, iYear, iAge) 
            S_Old = S_Source(iSplit, iAge)
            S_Smooth(iSplit, iAge) = (1 - S_Smooth_Factor) * S_Old + S_Smooth_Factor * S_Smooth(iSplit, iAge)
         enddo
      enddo
      
      ! Calculate residuals
      TotResSqr = 0
      do iYear = 0, Num_Years - 1
         do iSplit = size(YearGroups) - 1, 0, -1 
            if (iYear >= YearGroups(iSPlit)) exit
         enddo
         do iAge = 0, Num_TrueAges - 1
            N0 = Opt%PopModel%Nay(iYear, iAge)
            F_Sep = S_Smooth(iSplit, iAge) * Fy_Smooth(iYear);
            M = PValue(Opt, cNatMort, 0, iYear, iAge) 
            Z = F_Sep + M 
            C_Old = GetObsDataValue(Opt, odCaNum, 0, iYear, iAge)
            ! C = F * N0 * (1 - exp(-Z)) / Z
            C_New = F_Sep * N0 * (1 - exp(-Z)) / Z
            Opt%PopModel%InternalCatchISVPA(iYear, iAge) = C_New
            Res = 0
            if ((C_Old > 0) .and. (C_New > 0)) then
               Res = Log(C_Old / C_New)
            endif
            ResSqr = Res * Res
            TotResSqr = TotResSqr + ResSqr
            if(Finished) then
               ResArr(iYear, iAge) = Res
               ResSqrArr(iYear, iAge) = ResSqr
            endif
         enddo
      enddo
      
      Convergence = max(max((TotResSqr - LastRes(0))**2, (TotResSqr - LastRes(1))**2), (TotResSqr - LastRes(2))**2) 
      if (Convergence < 1e-5) exit
      LastRes(0) = LastRes(1)
      LastRes(1) = LastRes(2)
      LastRes(2) = TotResSqr
      
   
      ! Replace source selection and year factor with smoothed 
      do iSplit = 0, size(YearGroups) - 1
         do iAge = 0, Num_TrueAges - 1
            S_Source(iSplit, iAge) = S_Smooth(iSplit, iAge)
         enddo
      enddo
      do iYear = 0, Num_Years - 1
         if (iYear < Num_Years - 1) then
            Fy_Source(iYear) = Fy_Smooth(iYear)
         endif
      enddo
   enddo
   

   ! Replace original parameters with the converged source
   !do iSplit = 0, size(YearGroups) - 1
   !   iYear = YearGroups(iSplit)
   !   do iAge = 0, Num_TrueAges - 1
   !      pParamElm => GetParamElm(Opt, cSelection, 0, iYear, iAge, .false.)
   !      call SetPValueEx(pParamElm, S_Source(iSplit, iAge))  
   !   enddo
   !enddo
   !do iYear = 0, Num_Years - 1
   !  if (iYear < Num_Years - 1) then
   !      pParamElm => GetParamElm(Opt, cFYearFactor, 0, iYear, 0, .false.)
   !      call SetPValueEx(pParamElm, Fy_Source(iYear))
   !   endif
   !enddo
   ! Do not normalize selection and f year factor here, then the tuning will go ape
   
   ! Deliver N and F according to choosen effort control 
   ! If effort control, F = Fy * S overall and N0 = N1 * exp(Z)
   if (Opt%IndexInput%ISVPAControlType == ISVPA_EffortControl) then
      do iYear = Num_Years - 1 - 1, 0, -1
         do iSplit = size(YearGroups) - 1, 0, -1 
            if (iYear >= YearGroups(iSPlit)) exit
         enddo
         do iAge = Num_TrueAges - 1 - 1, 0, -1
            Fy = Fy_Smooth(iYear)
            S = S_Smooth(iSplit, iAge)  
            F = Fy * S
            N1 = Opt%PopModel%Nay(iYear + 1, iAge + 1)
            M = PValue(Opt, cNatMort, 0, iYear, iAge) 
            Z = F + M 
            N0 = N1 * exp(Z)
            call SetFay(Opt, iYear, iAge, F)
            call SetNay(Opt, iYear, iAge, N0)
         enddo
      enddo
   endif
   
   ! Stream out internal isvpa residuals
   !if(Finished) then
   !   open(cFileUnit, file = GetWorkFileName(Opt, "ISVPA_Residuals.txt"), action = 'write', status = 'replace')
   !   write(cFileUnit, "(A)") "ISVPA internal residuals = ln(obs/mod)"
   !   call print_matrix(ResArr, out_unit=cFileUnit, col_start = Opt%IndexInput%FirstYear, &
   !      row_start = Opt%IndexInput%FirstAge, &
   !      col_width = 18, &
   !      num_dec = 4 &
   !      )
   !   write(*, *) 
   !   write(cFileUnit, "(A)") "ISVPA internal residuals squared = ln(obs/mod)^2"
   !   call print_matrix(ResSqrArr, out_unit=cFileUnit, col_start = Opt%IndexInput%FirstYear, &
   !      row_start = Opt%IndexInput%FirstAge, &
   !      col_width = 18, &
   !      num_dec = 4 &
   !      )
   !   close(cFileUnit)      
   !endif
   
   deallocate(S_Source)
   deallocate(Fy_Source)
   deallocate(S_Smooth)
   deallocate(Fy_Smooth)
   deallocate(YearGroups)
end subroutine CalculateISVPAPopulation


subroutine OnCalculateFinished(Opt)
!   Purpose: Calculates the age plus group
!   Note:
!   History: When        Who        What
!            25-04-07    Skalevik   Created
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iYear, iStartYear, iAge;
   double precision :: C, N, F, M, Z, Na, Fa, Ma, Za, Np, Fp, Mp, Zp, Fy1, Fy2
   logical :: ByCatch
   if (NumAges(Opt%IndexInput) == 1) then
      return
   endif
   
   ! Year projection calculation 
   ! N and F are projected from a model spesific start year into the end of the model year (survivor stadium)
   ! ISVPA and VPA have startpoint set to last catch year + 1, separable have model year + 1 as startpoint.
   ! Enpoint is model year + 1 (End of model year = survivor stage)
   ! The projection goes up to the oldest true age, not including the plus group.
   iStartYear = NumBackwardProjectedYears(Opt%IndexInput)
   do iYear = iStartYear, NumYears(Opt%IndexInput)
      do iAge = 0, NumTrueAges(Opt%IndexInput) - 1
         select case(iAge)
            case (0)
               if (iYear < NumYears(Opt%IndexInput)) then 
                  N = PValue(Opt, cNValues, 0, iYear, 0)
                  ! If N value at parameter list is 0, it is most likely that it is not
                  ! tunable, hence we prefer -1 as a code to show that it is not taken
                  ! into account and the observed data are not playing a role here.
                  if(N == 0.0d0) N = -1.0d0
               else
                  N = -1.0d0
               endif
            case default
               N = Opt%PopModel%Nay(iYear - 1, iAge - 1)
               if(N > 0.0d0) then
                  F = Opt%PopModel%Fay(iYear - 1, iAge - 1)
                  if (F == -1) F = 0
                  M = PValue(Opt, cNatMort, 0, iYear - 1, iAge - 1) 
                  Z = F + M
                  N = N * exp(-Z)
               endif
         end select
         call SetNay(Opt, iYear, iAge, N)
         F = Opt%PopModel%Fay(iYear - 1, iAge)
         if (iYear < NumYears(Opt%IndexInput)) then
            if (F == -1) F = 0
            Fy1 = PValue(Opt, cFYearFactor, 0, iYear, 0)
            Fy2 = PValue(Opt, cFYearFactor, 0, iYear - 1, 0)
            if (Fy2 /= 0.0d0) then
               F = F * Fy1 / Fy2
            else
               F = -1.0d0
            endif
         endif
         call SetFay(Opt, iYear, iAge, F)
      enddo
   enddo

   ! Plus group calculation  
   iAge = NumAges(Opt%IndexInput) - 1
   ! 1. In plus group, set F as in the oldest true age
   do iYear = 0, NumYears(Opt%IndexInput)
      call SetFay(Opt, iYear, iAge, Opt%PopModel%Fay(iYear, iAge - 1))
   enddo

   do iYear = 0, NumYears(Opt%IndexInput)
      ByCatch = (iYear <= NumBackwardProjectedYears(Opt%IndexInput) - 1) .and. &
         ((0 == iYear) .or. (pgDirectFromCatch == Opt%IndexInput%PlusGroupMethod))
      if (ByCatch) then
         ! DIrect from catch
         C = GetObsDataValue(Opt, odCaNum, 0, iYear, iAge)
         M = PValue(Opt, cNatMort, 0, iYear, iAge) 
         F = Opt%PopModel%Fay(iYear, iAge)
         Z = F + M
         if ((F > 0) .and. (C > 0)) then
            N = C * Z / F / (1 - exp(-Z))
         else
            N = 0
         endif
      else 
         ! Dynamic pool (sum based on oldest age and plusgroup previous year).
         Na = Opt%PopModel%Nay(iYear - 1, iAge - 1)
         Fa = Opt%PopModel%Fay(iYear - 1, iAge - 1)
         if (Fa == -1) Fa = 0
         Ma = PValue(Opt, cNatMort, 0, iYear - 1, iAge - 1) 
         Za = Fa + Ma
         Np = Opt%PopModel%Nay(iYear - 1, iAge)
         Fp = Opt%PopModel%Fay(iYear - 1, iAge)
         if (Fp == -1) Fp = 0
         Mp = PValue(Opt, cNatMort, 0, iYear - 1, iAge) 
         Zp = Fp + Mp
         N = Na * exp(-Za) + Np * exp(-Zp)
      endif
      call SetNay(Opt, iYear, iAge, N)
   enddo

end subroutine OnCalculateFinished


end module Population
